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ary systems  which  are  driven  with  a  unit  step  as  the  test 
input.   Discrete  realization  of  the  continuous  systems  is 
obtained  using  the  measured  output-samples,  to  a  step  in- 
put, directly  in  the  algorithm. 
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I.       REVIEW    IF   HP's    ALGORITHM 

A.  INTRODUCTION 

Design  of  a  control  system  which  will  perform  accept- 
ably with  respect  to  some  criterion  and  satisfy  the  pos- 
sible constraints  over  its  operating  range  depends  on 
knowledge  of  the  process  dynamics.   The  dynamics  might  be 
linear,  nonlinear,  time  variant,  stationary  or  might  be  a 
function  of  environment.   In  any  case  the  dynamics  of  the 
process  must  be  formulated  by  a  set  of  differential  equa- 
tions (if  it  is  a  continuously  operating  system)  or  by  a 
set  of  difference  equations  (if  it  is  a  discrete  time  sys- 
tem) o      An  effective  design  and/or  analysis  can  proceed 
after  this  formulation. 

These  equations  might  be  completely  known  or  might  be 
written  down  from  the  parameter  values  supplied  by  the 
manufacturer  using  the  laws  of  physics.   Sometimes  only 
partial  information  or  no  information  is  supplied  about 
the  process ;   yet  somehow  the  dynamics  must  be  formulated. 
This  leads  to  the  solution  of  an  identification  problem. 

B.  BACKGROUND 

A  number  of  algorithms  can  be  found  [1,2,3,5]   for 
the  solution  of  the  identification  problem.   These  range 
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Numbers  in  brackets  indicate  the  references  at  the 
end. 


from  a  partial  identification,  like  the  damping  factor  of  a 
dominant  complex  pair  of  poles,  to  a  complete  identifica- 
tion of  the  systems  dynamics  using  frequency  and  time  do- 
main techniques. 

The  identification  problem  can  be  defined  in  general 
as:   to  find  an  internal  description  of  the  system  (a  set 
of  differential  or  difference  equations)  from  the  given 
external  description  (input-output  relation) . 

Input-output  relations  may  be  defined  in  the  time  do- 
main (impulse  response  for  continuous  systems  and  pulse 
response  for  discrete  systems)  or  in  the  frequency  domain 
(transfer  function).   In  the  identification  problem  this 
relation  is  given  in  general  as  experimental  data. 

Most  of  the  algorithms  are  applicable  to  linear,  sta- 
tionary systems  only;  and  only  this  class  of  systems  will 
be  considered  in  this  study. 

The  following  definitions  and  notations  will  be  used. 

A  continuous,  linear,  stationary  system  will  be  re- 
presented by 

x  =  Ax  +  Bu 

-^_         o/         ~ 

y  =  Cx 
where, 

A   is  a   nxn   system  matrix 

B   is  a   nxm   distribution  matrix 


(D2 


2  d 

x  denotes  -rr-  x 
~  dt  ~ 


C   is  a   pxn   observation  matrix 

x   =  x(t)   is  a   nxl   state  vector 

u   =   u(t)   is  a   mxl   input  vector 

Y  =  y(t)   is  a  pxl   output  vector 
and  a  discrete,  linear,  stationary  system  will  be  repre- 
sented by 

x(k+l)  =  A   x(k)  +  B   u(k) 

(2) 
y(k)  =  CD  x(k) 

with  dimensions  same  as  (1) 

A  continuous  system  is  said  to  be  controllable  [2]  if 
the  matrix 

\b\   ABl  A2B  I  ....  !  An_1B]  (3)3 

has  rank  n;   and  is  said  to  be  observable  if  the  matrix 

[cT!  ATCT  !  (AT)2CT!  ••-•  i  (AV-Vl  (4)4 

has  rank  n.   Same  definitions  hold  for  a  discrete  system  by 
replacing  A,  B,  C  with  A  ,  B_,  C   respectively. 

A  system  may  always  be  partitioned  into  four  possible 
subsystems  [2,3]  as  shown  in  Fig.  1. 

Part  I:     Controllable  and  observable 
Part  II:     Uncontrollable  but  observable 
Part  III:    Controllable  but  unobservable 
Part  IV:     Uncontrollable  and  unobservable. 


3 
Vertical  dotted  lines  indicate  partitioning. 

4  T 
C   indicates  transpose  of  C. 


u 


u 


-=*> 


II 


III 


IV 


Fig.  1 

As  is  seen  from  Fig.  1  only  Part  I  gives  the  input- 
output  relation.   Although  Part  II  seems  to  be  contribut- 
ing to  the  output,  if  no  energy  is  stored  at  this  part 
(zero  initial  condition)  any  contribution  will  be  the  re- 
sult of  a  noise  input  only. 

If  the  system  is  of  the  order  n  and  the  subsystems 


are  of  order  n  ,  n   ,  n    ,  n    respectively  then 


n  =  n_  +  n__  +  n___  +  n_.T7 
I     II     III    IV 


(5) 


Identification  (alternatively  called  realization)  of 
the  system  in  Fig.  1  will  give  a  system  of  order  n  . 

The  above  statements  are  given  as  the  following 
theorem  in  [2],  "Knowledge  of  the  impulse  response  identi- 
fies the  completely  controllable  and  completely  observable 
part,  and  this  part  alone,  of  the  dynamical  system  which 
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generated  it.   This  part  is  itself  a  dynamical  system  and 
has  the  smallest  dimension  among  all  realizations.   More- 
over, this  part  is  identified  by  its  impulse  response 
uniquely  up  to  algebraic  equivalence". 

Therefore  complete  identification  of  a  dynamical  sys- 
tem necessitates  complete  controllability  and  observability 
of  the  system.   Part  I  must  constitute  all  the  system  and 
parts  II,  III  and  IV  must  be  missing. 

In  this  study  only  completely  controllable  and  observ- 
able systems  will  be  considered. 

C.    Ho's  ALGORITHM  [l,4] 

First  the  algorithm  will  be  given  for  a  discrete  sys- 
tem as  represented  in  (2) . 

Let  h(k) . .  be  the  present  response  at  output  i  to  a 
unit  pulse  at  input  j  applied  k  time  units  ago.   Let  H,  be 
the  matrix 


\ 


h(k)l;L   h(k)12 


h(k)21   h(*)22 


h(k) 


pi 


h(k) 


im 


....  h(k) 


pm 


(6) 


Let  HO.  be  the  block  matrix 


HO,  = 


Hl     H2 


H2     H3 


H 


H 


,H 


21-1 


(7) 


Let  p.  be  the  rank  of  HO..   Build  the  block  matrices 

HO.  &   =  1,2,3,  ...  and  each  time  find  the  rank  p.,  until 

p .    equal  p.    . ,  let  this  value  of  p.   be  n  and  I   be  q;   then 

find  the  elementary  transformation  matrices  P  and  Q  such 

that 


P(HOq)Q  = 


1     0 


L° 


n  l's 


0 


(8) 
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Save  P  and  Q  matrices.   Let  HOT  be  the  block  matrix 

q 


HOT   = 

q 


H2     H3 


H3     H4 


11 


q+1 


II 


q+1 


H 


2q 


(9) 


Then  a  state  variable  description  of  the  system  is 

found  by  letting 

A^  =  the   nxn   northwest  corner  of   P(HOT  )Q 
D  q 

B^  =  the   nxm  northwest  corner  of   P(HO  )       (10) 
D  q' 

C^  =  the   pxn   northwest  corner  of   (HO  )Q 
D  q 

The  following  results  can  be  stated 

1.  p  =  p  , ,  =  n  is  the  minimum  dimension  for  the 

q    q+1 

realization. 

2.  By  choosing  P  and  Q  suitably  all  possible  realiza- 
tions may  be  obtained. 

3.  Any  two  minimal  realizations  are  isomorphic,  that 
is  there  exists  a  nonsingular  matrix  W  such  that 


An   =  W  A^   W 
D2       Dx 


-1 


B    =  W  B 
u2  ul 


(11) 


c   =  c   w 

D      D 
U2  Ul 


-1 
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and  W  is  given  explicitly  as 

W  =  (V2V2T)    V2V1T  (12) 


where  V   r  =  1,2   is  given  by 


V   = 
r 


'T  !  AT  CT  ;  (AT  )V   !  -..  ''(AT  )q_1CT 
r  '   r   i    r    r  !        r      r 


For  continuous  systems  the  algorithm  follows  the  same 
pattern  except  in  (6)  each  h(k) . .  must  be  replaced  by  the 
impulse  response  and  its  successive  derivatives  evaluated 
at   t  =  0,  that  is,  if  g(t) . .  is  defined  as  the  impulse 
response  between  input  j  and  output  i  then 

h(l)  .  .  =  g(0)  .  . 

11(2)..  =  g(0)±.  (14) 

h(3)  .  .  =  g(0)  .  . 

±3  13 


or  equally  g(t) . .  must  be  expanded  in  a  power  series  for 
all  t  as 

g(t)  .  .  =  £  ak  tk_1/(k-l):  (15) 

k=l 

and  each  h(k)  .  .  must  be  replaced  by  a,  .   These  necessitate 
2.J  K 

that  g(t) . .  must  be  a  real  analytic  function. 

If  the  input-output  relation  is  given  in  the  trans- 
form domain  (H(z) . .  for  discrete  systems  and  H(s) . .  for 
continuous  systems)  it  must  be  representable  as  a  power 
series  in  the  form  of 
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II 


(z)  .  .  =   )  b   z     (s  for  cont.  systems) 


(16) 


k=l 


and  each  h(k)  .  .  in  (6)  must  be  replaced  by  b,  ,  after  that 
the  algorithm  is  developed  in  similar  fashion. 

In  any  case  H,  ,  k  =  1,2,  ...  constitutes  a  sequence  of 
constant  matrices.   In  a  broad  sense  the  identification 
problem  may  be  restated  as  [l]:   Given  a  sequence  of  pxm 
constant  matrices,  H,  ,  k  =  1,2,  ...  find  the  triple 
A  ,  B  ,  C   (or  A,  B,  C)  of  constant  matrices  such  that 

Hk  =  CDADk"lBD  '   *'-!.*.  —  <"> 

to  show  the  mechanics  of  the  algorithm  the  following  ex- 
amples are  presented. 

1.   Example  1    Single -Input  j  Single-Output 

Discrete  System  (m  =  1,  p  =  1) 
Let  the  H,   k  =  1,2,  ...  be 

H,  =  1     H2  =    °     H3  =  "2 

H„  =  6     Hc  =  -14     Hc  =  30  ... 
4  5  6 

Start  building  HO.,  I   ~  1,2,  ... 

hoi  =  Chi  J  -  [ 1 1  •  pi  ■ x 


H02  = 


Hl    H2 


H2    H3 


0 


;   P2  =  2 
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HO. 


10-2 

0-2      6 

-2      6     -14 


P3  -  2 


therefore 

n  =  2   and   q  =  2 


HOT. 


H2     H3 


H3     H4 


0     -2 


-2 


A  set  of  P  and  Q  matrices  which  satisfy  (8)  is 


,   Q- 


0     -h 


P1(HOT2)Q1 


-2     -3 


=   A, 


P^HO^   = 


0     -2 


y      B 


D, 


(H02)Q1 


;Dl  ■  C1    °] 


^>    n 

1 
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therefore    the    state   equations    are 


x(k+l)       = 


xM   + 


-2  -3    ! 


u(k) 


0 


(k)  1  o]  x(k) 


this  realization  must  satisfy  (17) 


«k  =  V°i1Bdi  '  k  =  X'2' 


k-1 

A     can  be  calculated  by  the  Cayley-Hamilton  technique 

Dl 


Ak  -1  =  a  i  +  a.  a 
Dl     °     1   Di 


OL 


a. 


-2an    a  -3a. 
1     o   1 


a 


a. 


"k 


-£     0] 


2a,    a  -3a., 
l     o   l 


=  a 


where  a   and  a,  can  be  found  as  follows: 
o       1 

Eigenvalues  of  A    are 

Dl 


X      =   -1,       X      =    -2    ;    therefore 
(X^*-1  =   (-l)k_1  =  aQ  +  a1X1  =  aQ  -  ax 
(X2)k_1  =    (-2)k_1  =  aQ  +  axX2  =  aQ  -20i± 
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Substituting  for  k  =  1,2,  .. 


1-1        1-1 
H1  =  2(-l)1  x  -  (-2)1  x  =  1 


2-1        2-1 
H2  =  2(-l)^  L    -  (-2T    =  0 


3-1         3-1 
H3  =  2(-l)J  L    -    (-2)J    =  -2 


which  are  exactly  the  same  as  the  starting  H,  ,  k  =  1,2, 
A  second  set  of  P  and  Q  matrices  is 


0 


,      Q. 


0 


which  give  the  realization 


x(k+l)  = 


0     -2 


-3 


x(k)  + 


u(k) 


(k)  =  [l      o]x(k) 


this  realization  also  satisfies  (17) . 

For  the  isomorphism  of  two  realizations,  W  is  calculated  as 


w  =  (v2v2T)  v2VlT 
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By    (13)  : 


V 


=    \CT      '  AT    CT      I    = 
±i         11 


0 


1  0 


0  1 


v^  = 


-2 


1  0 


0         -2 


W      = 


1  0 


o       -h 


;  W 


-1 


1  0 


0         -2 


A^      =   WA„  W 
D2  Dl 


-1 


1  0 


0  -h 


0  1 


•2         -3 


1  0 


0         -2 


0         -2 


1         -3 


D2  D1 


1  0 


0         -h 


0 


CD2    -   CDlW_1       =    t  °] 


1  0 


0         -2 


-"&      o] 
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The  two  realizations  have  the  following  signal  flow 
graphs : 


U(z) 


•  Y(z) 


Y(z) 


U(z) 


Fig.  2 


By  Mason's    Gain   Formula  both  have    the    transfer   function 


H(z)    -         (z+3) 
tiKZ}    ~    (z+1)  (z+2) 


By  a  long  division 

H(z)  =  z"1  +  0z~2  -  2z~3  +  6z-4 

-k 
the  coefficients  of  z    ,  k  =  1,2,  . , 

H,  ,  k  =  1/2,  ... 


-  14z 


-5 


are  again  equal  to 


18 


2.   Example  2    Two— Inputs  j    Three-Outputs 

Continuous  System  (m  =  2,  p  =  3) 
Le t  the  H,  ,  k  =  1,2,  ...  be 


H, 


H„   = 


1 

0 

0 

1 

1 

1 

1 

2 

'  H2  " 

-2 

-4 

•  H3  ~ 

4 

8 

2 

1 

-3 

-4 

8 

14 

3 

-5 

7 

13 

, 

8 

-16 

-  H5  = 

16 

32 

18 

-34 

38 

74 

. 

Then 


HO.  = 


H02  = 


1 

0 

1 

2 

'    Pl 

=  2 

2 

i 

1_ 

1 

0 

0 

1 

1 

2 

-2 

-4 

2 

1 

-3 

-4 

0 

1 

1 

1 

-2 

-4 

4 

8 

-3 

-4 

8 

14 

P->  =  3 
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;  i 

1 

,  0 

0 

1 

1 

1 

1 

2 

-2 

-4 

4 

8 

2 

1 

-3 

-4 

8 

14 

0 

1 

1 

1 

-3 

-5 

HO3   = 

-2 

-4 

4 

8 

-8 

-16 

p3  = 

=  3 

-3 

-4 

8 

14 

-18 

-34 

•' 

r 

1 

-3 

-5 

7 

13 

1 

1 

4 

8 

-8 

-16 

16 

32 

8 

L_ 

14 

-18 

-34 

38 

74 

therefore 

n   =    3         and        q  =   2 


HOT2    = 


0 

1 

1 

2 

-4 

4 

3 

-4 

8 

1 

1 

-3 

4 

8 

-8 

8 

14 

-18 

ll 

8 

14 

-5 

-16 

■34 
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A  set  of  P  and  Q  matrices  is 


Pl  = 


Ql  = 


1 

0 

-2 

2 

-2 

-1 

1 
0 

0 
0 


0 
0 
0 

1 

0 
0 

0 

1 

0 
0 


0 

0 

0 

0 

0 

1 

0 

0 

1 

-1 

0 

0 

2 

0 

1 

-1 

2 

2 

1 

0 

1 

-1 

-1 

1 

3 

-4 

2 

3 

5 

-7 

3 

4 

which  give  the  realization 

0      10 

111 

-4     -7     -4 


x   + 


0 

1 

0 


u 


y  = 


10  0 
12  1 
2      11 


x 
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this  has  the  impulse  response  matrix 


G(t) 


3/2  -  €-fc  +  §  r2t 


-2t 


-t      -2t 
2-3e    +  2€ 


2e 


-2t 


3/2  -2€  fc  +  5/2e  2t     2-6€  fc  +  5c  2t 


Expanding  the  exponentials  in  power  series  and  expressing 
the  result  as  a  matrix  polynomial  the  following  expression 
can  be  obtained. 


G(t) 


1  0 

1  2 

2  1 

-3  -5 

-8  -16 

-18  -34 


0 

1 

2 

-4 

t    + 

3 

-4 

1 

1 

4 

8 

8 

14 

t2/2 


X31 


7  13 
16  32 
38      74 


4 

b/4:  + 


Each  element  of  the  coefficient  matrices  corresponds  to  the 
a,  '  s  of  (15)  and  the  matrices  themselves  are  exactly  H,  's. 
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A  second  set  of  P  and  Q  matrices  is 


1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

1 

0 

1 

0 

0 

0 

2 

1 

-2 

0 

1 

-1 

2 

0 

2 

2 

1 

0 

1 

0 

1 

-1 

-1 

1 

-2 

-3 

3 

-4 

2 

3 

-2 

3 

Q 

2    ~ 

-5 

-5 

5 

-7 

3 

3 

-3 

4 

which  give  the  realization 


x   = 


0 

1 

0 

0 

0 

1 

0 

-2 

-3 

X    + 


1 

0 
1 


0 

1 
1 


u 


y  = 


10  0 
Oil 
10     1 


x 
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It  can  be  shown  as  before  that  this  realization  also 
satisfies  (17) 

To  show  the  isomorphism  of  two  realizations  W  is 
calculated  as 


i "  |_Ci  i Ai  Ci  _ 


V.  = 


1  1  i  2  0-2-3 
0  2  1  1-4-4 
0     1     1     0-2-3 


V2  = 


10  1 
0  10 
0     11 


0  0     0 

1  -2    -1 
0    -2    -3 


41 

8 

-    7 

T    ~1             1 

(v2v2    )        =  ft 

8 

29 

-16 

-7 

-16 

14 

v2vx- 


3      11 

8     15      8 

16     23     15 
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m      -1  m 

w       =   (v2v2T)     v2v±T 


10  0 

0  10 

111 


w 


-1 


10  0 

0  10 

-1      -1        1 


1 

0 

0 

0 

1 

1 

a 

1 

0 

0 

A2 

=  WA..W    1    = 

0 

1 

0 

1 

1        1 

0           1 

1 

1 

1 

1 

-4 

-7         -4 

[_-l         -1 

1 

0 

1 

0 

0 
0 

0 
-2 

1 

-3_ 

1 

0 

0 

1 

0 

1           0 

B2 

=   WB 

0 

1 

0 

0 

1 

= 

0           1 

1 

1 

1_ 

_0 

0_ 

_  1        1_ 

, 

__ 

— i     i — 

1 

0 

0 

1 

0           0 

1           0 

0 

C2 

=    CJN           = 

1 

2 

1 

0 

1           0 

= 

0           1 

1 

2 

1 

1 

-1 

-1 

L 

1 

0 

1 
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The  W  matrix  also  relates  the  two  sets  of  P  and  Q 
matrices  as 


w 

0 

V1 

0 

P2    " 

pl  ' 

'       Q2    = 

=  Q1 

0 

I 

0 

I 

— 

— 

For  this  specific  example 


1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

P2    " 

1 

1 

1 

0 

0 

0 

' 

. 

0 

0 

0 

1 

0 

0 

. 

0 

1 

0 

0 

0 

1 

0 

0     • 

0 

0 

0 

0 

1 

1 

0 

0 

0 

0 

1 

0 

0 

Q2 

=      Q1 

-1 

-1 

1 

0 

0 

0 

0 

1 

D.    ANOTHER -ALGORITHM 

The  algorithm  given  in  [5]  for  discrete  systems 
(equally  applicable  to  continuous  systems)  with  scalar 
measurement  (single— output)  shows  a  similarity  with  Ho ' s 
Algorithm. 
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For  the  algorithm,  a  free  dynamical  system  (no  forc- 
ing function)  represented  as  in  (18)  is  considered 


x(k+l)  =  A   x(k) 
~         jj  — 


y(k)    =  CD  x(k) 


(18) 


where 

A    is'  a   nxn  matrix 

C    is  a   lxn   row  vector 

x    is  a   nxl   state  vector 

y    is  the  scalar  output 
and  x(0)  is  known.   Then  HO  and  HOT.  matrices  are  construe- 
ted  as  in  Ho ' s  Algorithm  for  single-input^ single-output 
systems,  from  the  measured  output— sequence.   An  algebraic 
equivalent  of  the  system  is  realized  by  letting 

A '  =  (HOT  ) (HO  )   ■ 
D        q     q 

(19) 
C^  =  [1   0   0] 

Note  that  for  single-input ^single-output  systems  of  order 
n,  q  is  equal  to  n,  and  p  also  is  equal  to  n.   Then  this 
algorithm  becomes  a  special  case  of  Ho ' s  Algorithm  by  let- 
ting 

P  =  I 

_1  (20) 

Q  =  (H0q) 

Because  a  nonsingular  matrix  can  always  be  reduced  to  an 
identity  matrix  by  elementary  transformations  on  its 
columns  (or  rows)  only. 
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v,    For  complete  identification  of  the  system,  complete 
controllability , and  observability  was  a  necessary  and  suf- 
ficient condition  for  systems  with  forcing.   Here,  as  B 
is  absent,  controllability  cannot  be  checked;   instead  a 
new  condition  called  n-identification  is  imposed. 

A  system  is  said  to  be  n-identifiable  if  the  matrix 
[X(0)  I  ADx(0)  J  AD2x(0)  J  ...  J  ADn_1x(0)]         (21) 

has  rank  n.   Similarity  of  (3)  and  (21)  is  obvious  for 
single-input  systems. 

Therefore  it  is  possible  to  determine  A   and  C    if 
the  system  is  observable  and  n-identifiable.   These  two 
conditions  are  summed  together  and  called- 1-identifiability 
[5]. 

In  fact  another  realization  is  possible  by  letting 

A^  =  (H0q)_1  (HOTq) 

(22) 

C '  =  [first  row  of  HO  ] 
D  q 

This  is  also  a  special  case  of  Ho's  Algorithm;   simply 

P  =  (HO)'1 

q  (23) 

Q   =   1 

The  following  relation 

AD  =  <AD)T  <24> 
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holds  for  these  two  extreme  cases.   As  HO  and  HOT  are  sym- 
metric matrices  for  single-input ,single-output  systems  it 
can  be  written 


(a')T  =  [(HOTq)(HOqr  f 

"    [(HV"1]T[H0Tq]T  <25> 

=       (H0q) -1     (HOTq)    =  A^ 

Example  1  is  applicable  to  these  cases  just  simply 
letting 

*<°)   =  BD 
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II.   APPLICATION  OF  Ho ' s  ALGORITHM 

A.    IMPLEMENTATION  OF  Ho ' s  ALGORITHM 

Ho ' s  Algorithm  perhaps  is  the  one  which  requires  min- 
imum computation  among  the  existing  identification  schemes. 
Because  of  the  simplicity  of  the  computations  it  can  be 
used  in  many  applications. 

The  algorithm  is  implemented  with  a  digital  computer 
program  to  see  how  it  can  be  used  in  an  off  line  fashion 
for  identification.   The  results  obtained  can  be  extended 

for  use  in  other  applications. 

5 
1.   Digital  Computer  Programs 

The  program  consists  of  one  main  program  and  four 
subroutines.   The  main  program  reads  the  data,  builds  the 
HO  and  HOT  matrices  and  controls  the  subroutines.   For  dig- 
ital simulation  of  systems,  data  can  be  generated  within 
the  main  program  by  replacing  the  data  reading  loop  with  a 
generating  function.   The  algorithm  does  not  stop  when  two 
successive  ranks  of  HO  are  equal,  but  uses  all  available 
data  for  additional  checking. 

Subroutine  SDWFD2  finds  the  successive  derivatives 
evaluated  at  zero  time  for  continuous  systems;   for  dis- 
crete systems  it  is  by-passed.   FLAG  controls  this  opera- 
tion. 


5 
See  pages  56-61 
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Four  Newton  Forward  Differences  are  used  to  evaluate  each 
derivative.   Subroutine  RANK  finds  the  rank  of  HO  each 
time  it  is  constructed.   Subroutine  PHOQ  finds  the  elemen- 
tary transformation  matrices  P  and  Q  which  put  HO  into  its 
normal  form.   Subroutine  MULT  makes  the  matrix  multiplica- 
tions required  by  the  algorithm. 

The  program  is  written  for  single  input-output  systems, 
but  can  be  altered  for  multiple  input-output  systems;   in 
fact  only  the  reading  data  and  building  HO  and  HOT  parts 
need  changes  in  the  main  program.   No  change  is  required 
in  the  subroutines. 

B.    APPLICATION  TO  CONTINUOUS  SYSTEMS 

For  discrete  systems  the  algorithm  is  straight  forward, 
the  only  limiting  factor  is  the  accuracy  of  the  measure- 
ments.  For  continuous  systems  however  some  difficulties 
arise.   First,  as  a  test  input  signal  an  impulse  is  needed. 
For  linear  systems  the  impulse  response  is  the  derivative 
of  the  response  to  a  step  input.   Therefore  a  step  input  can 
be  used  as  the  test  input.   Higher  order  inputs  can  also  be 
used  if  it  is  necessary  as  long  as  the  right  derivatives 
are  obtained.   The  program  is  written  to  accept  the  step 
response  of  the  system. 

A  second  problem  comes  from  the  numerical  evaluation 
of  successive  derivatives  from  the  measured  samples. 
Numerical  evaluation  of  a  derivative  is  an  approximation 
even  with  exact  measurements,  when  higher  derivatives  are 
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needed  accuracy  diminishes  greatly.   A  small  error  in  the 
input  data  deteriorates  the  results,  as  the  numerical  meth- 
ods magnify  these  errors  enormously.   The  sampling  period 
is  the  biggest  factor  in  this  magnification  as  it  appears 
in  the  denominator,  and  for  each  additional  higher  deriva- 
tive its  power  is  increased  by  one.   Magnification  of  error 
increases  as  sampling  period  decreases;   but  on  the  other 
hand  for  fast  changing  functions  a  small  sampling  period 
is  needed. 

The  examples  presented  in  the  following  tables  were 
simulated  digitally  for  unit  step-inputs  to  the  systems. 
Two  parameters,  DX  and  DD,  in  addition  to  the  sampling 
period,  seemed  to  have  great  effect  on  the  accuracy  of  the 
obtained  results.   A  suitable  choice  of  the  three  gave  re- 
pea ting  answers  when  the  dimension  of  the  HO  matrix  became 
larger  than  the  order  of  the  system.   DX  and  DD  are  posi- 
tive small  numbers;   when  finding  the  rank  of  HO,  any 
number  which  has  an  absolute  value  less  than  DX  is  equated 
to  zero  by  the  program.   When  finding  the  P  and  Q  matrices, 
which  put  HO  to  its  normal  form,  any  number  which  has  an 
absolute  value  less  than  DD  is  made  equal  to  zero  by  the 
program.   The  values  of  these  parameters  which  gave  a  good 
result  for  this  specific  example  are  included  in  the 
tables. 

In  example  3  a  system  with  two  real  poles  is  consider- 
ed.  The  poles  are  not  greatly  separated  and  a  gc od  real- 
ization became  possible  because  the  step  response  is  a 
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fairly  smooth  function,  and  in  the  evaluation  of  the  deri- 
vatives not  too  much  error  is  involved,  in  fact  the  first 
four  derivatives  are  almost  exact.   In  Example  4  a  system 
with  three  real  poles,  which  are  separated  by  a  significant 
amount,  is  considered.   In  part  A,  the  sampling  rate  is  not 
so  high.   Identification  of  the  pole   at  the  origin  is 
almost  exact,  the  pole  at  -1000  is  identified  fairly  well 
but  the  far  pole  is  in  error  by  a  large  amount.   In  part  B, 
the  sampling  rate  is  increased;   while  the  pole  at  -1000  is 
identified  more  accurately  and  the  far  pole  identification 
showed  improvement,  identification  of  the  pole  at  the  ori- 
gin began  to  deteriorate.   In  Example  5,  a  highly  damped 
complex  pair  of  poles  gave  a  very  accurate  identification 
with  a  high  sampling  rate;   as  there  is  no  spread  of  poles 
a  good  choice  of  sampling  rate  results  in  an  accurate  real- 
ization.  In  Example  6,  a  very  lightly  damped  complex  pole 
pair  is  considered.   Again  identification  is  very  accurate 
because  of  the  same  reasons  in  Example  5. 

In  the  above  examples  simulation  of  systems  was  carried 
out  digitally  in  double-precision.   In  general  a  high  samp- 
ling rate  with  a  high— precision  measurement  is  necessary 
which  makes  the  scheme  impractical  to  use  with  measurement 
of  analog  signals.   A  different  approach  is  necessary  which 
is  the  subject  of  the  next  section. 
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1 

EXAMPLE  3 

Transfer 
Function 
of  The 
System 

300 

(s+5) (s+20) 

Sampling  Per. 

0.005 

DX 

0.005 

DD 

0.000001 

Numerically 
Evaluated 
Successive 
Derivatives 

0.3376xl0"3     0.3000xl03    -0.7500xl04 

0.1575xl06     -0.3187xl07     0.6393xl08 

-0.1278xl010     0.2448X1011    0.2830xl013 

Realized 
Sys  tern 
Matrix 

-0.2132xl0~13       -O.lOOOxlO3 
0.1000x10          -0.2500xl02 

Realized 
Distribution 
Matrix  T 

0.1000x10           0.2541x10   u 

Realized 

Observation 

Matrix 

0.3376xl0~3         0.3000xl03 

Eigenvalues 
of  System 
Matrix 

-0.4999  x  101  +  j  0.0 
-0.1999  x  10   +  j  0.0 

34 


EXAMPLE  4a 

Transfer 
Function 
of  The 
System 

io12 

s(s+1000) (s+1000000) 

Sampling  Per. 

0.001 

DX 

1.0 

DD 

1 

1.0/1  largest  element  of  HO 

1  x  2.0 

1 

Numerically 
Evaluated 
Successive 
Derivatives 

0.430xl02     0.9148xl06 

0.8361xl012  -0.7993xl015 

21            24 
-0.7304x10     0.6981x10 

0.8746xl09 

1  p 
0.7641x10 

-0.6670xl027 

Realized 

System 

Matrix 

-0 . 5684xl0~13 -0 . 1309xl0~9 
O.lOOOxlO1    0.3725xl0~8 
0.1065xl0~13  O.lOOOxlO1 

0.2131xl0_1 
-0.1991xl07 
-0.3039xl04 

Realized 
Distribution 
Matrix  T 

O.lOOOxlO1    0.3552xl0-14 

-0.1387xl0"16 

Realized 

Observation 

Matrix 

0.430xl02     0.9148xl06 

-0.8746xl09 

Eigenvalues 
of  System 
Matrix 

-0.2083  x  104  +  j  0.0 

-0.956   x  103  +  j  0.0 

0.1070  x  10_7+  j  0.0 
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EXAMPLE  4b 


Transfer 
Function 
of  The 
System 


10 


12 


s(s+1000) (s+1000000) 


Sampling  Per, 


0.0001 


DX 


1.0 


DD 


1.0/1  Largest  element  of  HO  |  x2.0 


Numerically 
Evaluated 
Successive 
Derivatives 


-0.9625 


0.1000x10 


-0.9918x10 


0.8123xl012   0.2927xl016   -0.8084xl020 


0.1704xl025  -0.3552xl029 


0.7400x10 
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Realized 

System 

Matrix 


0.8633x10  12-0. 5544x10  8   -0.5187xl03 
O.lOOOxlO1   -0.1989xl0_11  -0.2083xl08 


0.6098x10  19  O.lOOOxlO1 


-0.2183x10' 


Realized 
Distribution 
Matrix  T 


O.lOOOxlO1   -0.5963x10  18  -0 . 2646xl0~22 


Realized 

Observation 

Matrix 


■0.9625 


0.1000x10 


-0.9918x10 


Eigenvalues 
of  System 
Matrix 


-0.2083  x  10  +  j  0.0 
-0.9999  x  103  +  j  0.0 
-0.2490  x  10_4+  j  0.0 
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EXAMPLE  5 

Transfer 
Function 
of  The 
System 

1000000 

(s+25000  +  J25000) (s+25000  -  J25000) 

Sampling  Per. 

0.000001 

DX 

0.000005 

DD 

1.0/1  Largest  element  of  HO  1  x  2.0 

Numerically 
Evaluated 
Successive 
Derivatives 

0.6764X101     0.8164xl06     -0.4927X1011 
0.1443xl06    -0.1056xl020    -0.1306xl025 
0.3408xl03    -0.2533xl037     0.2595xl044 

Realized 
Sys  tern 
Matrix 

0.9094xl0-11     -0.1250xl010 
O.lOOOxlO1       -0.5000xl05 

Realized 
Distribution 
Matrix  T 

O.lOOOxlO1        0.1058xl0-20 

Realized 

Observation 

Matrix 

0.6764X101        0.8164x10 

Eigenvalues 
of  System 
Matrix 

-0.2500  x  105  +  j  0.2500  x  105 
-0.2500  x  105  -  j  0.2500  x  105 
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EXAMPLE  6 


Transfer 
Function 
of  The 
Sys  tern 


1001 


(s  +  0.01  +  J500) (s  +  0.01  -  j500) 


Sampling  Per. 


0.0001 


DX 


0.005 


DD 


0.000001 


Numerically- 
Evaluated 
Successive 
Derivatives 


-0.2084x10 


-6 


■O.lOOlxlO4         -0.1986xl02 


-0.2502x10 


0.2890x10 


17 


0.1128x10 


-0.3539x10 


8 


22 


0.6234x10 


0.4010x10 


14 


27 


Realized 

System 

Matrix 


0.8198x10 


-13 


-0.2500x10 


0.1000x10" 


-0.1984x10 


-1 


Realized 
Dis  tribution 
Matrix  ™ 


0.1000x10' 


-0.2082x10 


-9 


Realized 

Observation 

Matrix 


-0.2084x10 


-6 


0.1001x10 


Eigenvalues 
of  System 
Matrix 


-0.9922  x  10  2  +  j  0.4999  x  103 
-0.9922  x  10~2  -  j  0.4999  x  103 
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C.    DISCRETE  REALIZATION  OF  CONTINUOUS  SYSTEMS 

As  is  seen  from  the  preceding  examples,  to  find  a  set 
of  differential  equations  for  continuous  systems  is  some- 
what a  trial  and  error  procedure,  as  one  must  come  up  with 
a  suitable  triple  of  DX,  DD  and  sampling  period,  if  no 
apriori  knowledge  is  available  for  the  system.   More  impor- 
tantly, very  accurate  measurements  are  needed  which  is  im- 
possible in  an  actual  experiment  with  hardware. 

The  main  error  source  in  the  above  procedure  is  the 
evaluation  of  derivatives  numerically.   If  the  measured 
data  can  be  used  directly,  without  evaluating  derivatives, 
inclusion  of  this  error  to  the  algorithm  can  be  eliminated. 
This  approach  leads  to  a  system  which  is  operating  con- 
tinuously but  observed  discretely  as  in  Fig.  3 


u(t) 


Continuous 
Sys  tern 


S   >*  y(kT) 


->.  y(t) 


Fig.  3 

where  T  is  the  sampling  period  and  u ( t)  is  a  step.   If  the 
system  has  the  dynamics  as  in  (26) 


x  =  Ax  +  Bu 
y  =  Cx 


(26) 
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Then 


x(t)  =1    {(Sl-A)  1   x(t  )  +  (SI-A)  X  Bu(s)}      (27) 


As  u(t)  is  a  step  (amplitude  V) 

x(t>  =X_1  {(sI-A)"1  x(t0)  +   (si-a)-^^ 

Defining 

0(t-tQ.)  =X~   {(sI-A)"1} 

r(t-to)  *£*  {(sI-a)"1b 


(28) 


(29) 


Substituting  in 

\    x(t)  =  0(t-t  )x(t  )+  T(t-t  )V  (30) 

^     '  *    o  ^     o  o 


t   =  kT 
o 

t   =  t   +  T 
o 


letting 

=  kT 

(31) 

x[(k+l)T]  =  0(T)x(kT)  +  T(T)V  (32) 

and 

y(kT)  =  Cx(kT)  (33) 

0(T)  and  F(T)  are  constant  matrices  as  the  sampling 
period  is  constant,  V  can  be  replaced  by  v(kT)  in  (32)  let- 
ting 

v(kT)  =  1   ,   k  =  0,1,2,  ...  (34) 

as  the  test  input  is  a  unit  step. 
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Dynamics  described  by  (32)  and  (33)  can  be  represented 
as  a  discrete  system  of  Fig.  4. 


v(kT) 


Discrete 
System 


y(kT) 

>» 


Fig.  4 

If  Fig.  3  is  driven  by  a  unit  step  and  Fig.  4  with  a 
unit  impulse  train  of  period  T,  they  both  will  give  the 
same  y(kT) .   Therefore  samples  taken  as  the  result  of  ex- 
periment on  Fig.  3  can  be  used  to  identify  the  system  in 
Fig.  4  as  if  it  is  driven  by  a  unit  impulse  train  and  the 
output  is  measured. 

Dynamics  described  in  (32)  and  (33)  is  an  exact  re- 
presentation of  the  continuous  system  of  Fig.  3  at  sampling 
instants  and  constitutes  a  discrete  realization  of  the  con- 
tinuous system.   Two  points  of  importance  must  be  mentioned 
here.   One  is  the  loss  of  controllability  and  observability 
due  to  the  sampling  of  continuous  systems  [2],  as  a  result 
a  realization  which  is  not  a  faithful  representation  of 
the  system  will  be  obtained.   If  a  linear-stationary— con- 
tinuous system  is  controllable  and  observable,  it  will  re- 
main so  when  sampling  is  introduced  if  and  only  if  [6] 

k2ff 


Re[s.]  =  Re[s.]  implies   Im[s.-s.]  ^ 


(35) 
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where 

s   ,  m  =  1,2,  ...n   eigenvalues  of  system 
m 

i  ,  j  =1,2,  ...  n 

k  =  positive  integer 

T  =  sampling  period 
If  no  apriori  knowledge  is  available  about  the  system,  at 
least  two  realizations  must  be  obtained  with  different 
sampling  rates  to  avoid  misrepresentation. 

The  second  point  comes  from  the  usage  of  a  unit  step 
as  the  test  input.   As  a  result  the  discrete  equivalent  in 
Fig.  4  is  driven  with  a  unit  impulse  train.   If  the  ori- 
ginal system  has  a  transfer  function 

H(z)  =  fjfj-  =  a1z_1  +  a2z"2  +  a3z"3  +  ...         (36) 
in  the  z-domain,  the  realization  will  be  for 

»'  <Z>  '=  U-l)D(z)  "  aiZ_1  +  a2Z"2  +  a3Z"3  +  •••   <"> 
a.  ,  i  =  1,2,  ...  are  measured  sample  values  to  a  unit 
step  input.   This  will  give  an  additional  eigenvalue  in- 
creasing the  order  of  the  system  by  one.   This  can  be 
avoided  by  a  simple  manipulation  of  the  data. 

H(z)  =  H'(z)(l-z_1)  =  a^z-1  +  (a^-a-^z-2  +  •••    (38) 
Therefore 

al  =  al 

a.     =  a'  -  aj  (39) 

a3  =  a3  "  a2 
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This  manipulation  of  the  data  has  an  additional  benefit,  if 
there  is  an  unknown  biasing  in  the  measurement  it  is  remov- 
ed automatically  except  from  the  first  sample.  The  program 
takes  care  of  this  conversion  of  the  data. 

In  the  following  examples  systems  are  simulated  digit- 
ally.  To  see  the  sensitivity  of  the  algorithm  to  numerical 
precision,  sample  values  were  rounded-off  to  four  signific- 
ant digits,  and  a  realization  without  manipulation  of  the 
data  is  included  in  some  examples  to  make  comparisons  pos- 
sible.  As  will  be  seen  the  choices  of  DX  and  DD  are  not 
critical  as  in  the  continuous  representation  of  systems. 

Example  7a  shows  the  identification  of  the  additional 
pole  at  the  origin  because  of  a  unit  step  test  signal.   In 
7b  this  pole  is  removed  with  the  manipulation  of  the  data, 
and  six  significant  digits  are  used  for  sample  values, 
realization  is  very  accurate.   In  7c  sample  values  are 
rounded  to  four  significant  digits,  which  results  in  a  de- 
terioration in  the  answer.   In  Example  8a  and  8b  seven 
significant  digits  are  used.   In  8a  an  additional  pole  is 
present  because  of  step  input,  in  8b  it  is  removed;   with 
the  same  significant  digits  8b  gave  a  more  accurate  answer 
as  the  system  order  decreased  by  one,  less  computation  is 
required  and  error  is  decreased.   In  8c  four  significant 
digits  are  used  with  some  deviation  of  the  answer  as  a  re- 
sult.  In  Example  9a  and  9b  a  complex  pole  pair  and  a  real 
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pole  are  considered,  with  eight  and  four  significant  digits 
of  sample  values  respectively.  Highly  accurate  results  are 
obtained  in  both  cases.  Example  10a  shows  the  loss  of  con- 
trollability and  observability  (as  the  sampling  is  intro- 
duced in  the  measurement,  only  observability  is  lost)  due 
to  the  sampling.  In  10b  with  a  different  sampling  rate  the 
system  is  identified  completely. 


■ 


■  • 
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EXAMPLE  7a 

Transfer 
Function 
of  The 
System 

160 

(s+1) (s+2) 

Sampling  Per. 

0.1 

DX 

0.001 

DD 

0.00001 

The  Type 
of  The  Data 

Not  manipulated  to  take  the  step  out 
Double-Precision 

Realized 

System 

Matrix 

0.3552xl0"14    0.0              0.7408 
O.lOOOxlO1      0.3552xl0"14    -0.2464X101 
0.6661xl0~15    O.lOOOxlO1      0.2723x10 

Realized 
Distribution 
Matrix  T 

1               -14              -15 
0.1000x10       0.3552x10       -0.2220x10 

Realized 

Observation 

Matrix 

0.7244         0.2628X101       0.5374X101 

Eigenvalues 
of  System 
Matrix 

-14 
0.8187  +  j  0.1230  x  10 

0.9048  -  j  0.7288  x  10"15 

0.1000  x  101  -  j  0.2778  x  10"15 

Corresponding 
s -Domain 
Eigenvalues 

-2.0 

-1.0 

0.0 
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EXAMPLE  7b 

Transfer 
Function 
of  The 
System 

160 

(s+1) (s+2) 

Sampling  Per. 

0.1 

DX 

0.001 

DD 

0.00001 

The  Type 
of  The  Data 

Manipulated  to  take  the  step  out 
Rounded  to  six  significant  digits 

Realized 
Sys  tern 
Matrix 

0.2220  x  10~15       -0.7408 
0.1000  x  101          0.1723  x  101 

- 

• 

1 

■  . 

Realized 
Distribution 
Matrix  T 

0.1000  x  101         -0.5551  x  10_16 

■ 

Realized 

Observation 

Matrix 

0.7244                0.1904  x  101 

Eigenvalues 

of  SVifpni 

0.8187  +  j  0.0 
0.9048  +  j  0.0 

Matrix 

Corresponding 
s -Domain 
Eigenvalues 

-2.0 
-1.0 
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EXAMPLE  7c 

Transfer 
Function 
of  The 
System 

160 

(s+1) (s 

+2) 

Sampling  Per. 

0.1 

DX 

0.001 

DD 

0.00001 

The  Type 
of  The  Data 

Manipulated  to 
Rounded  to  four 

take  the  step  out 
significant  digits 

Realized 

System 

Matrix 

0.6661  x  10~15 
0.1000  x  101 

-0.7386 
0.1722  x  101 

Realized 
Distribution 
Matrix  T 

0.1000  x  101 

-0.6938  x  10"16 

Realized 

Observation 

Matrix 

0.7245 

0.1904  x  101 

Eigenvalues 
of  System 
Matrix 

0.8069  +  j  0.0 
0.9154  +  j  0.0 

Corresponding 
s -Domain 
Eigenvalues 

-2.145 
-0.885 
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EXAMPLE  8a 

Transfer 
Function 
of  The 
System 

4000 

(s+1) (s+5) (s+10) 

Sampling  Per. 

0.1 

DX 

0.001 

DD 

0.00001 

The  Type 
of  The  Data 

Not  manipulated  to  take  the  step  out 
Rounded  to  seven  significant  digits 

Realized 
Sys  tern 

0. 426x10" 13  -0.5684xl0-13  -0  .  5684xl013  -0  .  2030 
O.lOOOxlO1  -0.5684xl0-13  -0.5684xl013  0.1310X101 

Matrix 

0.0         O.lOOOxlO1    0.5684xl013-0.2987xl01 
-0.3907xl0~13  0.0           O.lOOOxlO1   0.2880x10 

Realized 
Distribution 
Matrix  T 

O.lOOOxlO1   0.0           0.0          0.0 

Realized 

Observation 

Matrix 

0.4537       0.2542X101    0.6169X101   0.1077xl02 

Eigenvalues 
of  System 
Matrix 

0.9987  -   j  0.4198  x  10-14 
0.9073  -   j  0.2743  x  10~13 
0.3715  -   j  0.2742  x  10_13 
0.6028  +   j  0.5936  x  10~13 

Correspond- 

0.0 

ing 

s-Domain 

Eigenvalues 

-0.972 

-9.9 

-5.06 
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EXAMPLE  8b 


Transfer 
Function 
of  The 
System 


4000 


(s+1) (s+5) (s+10) 


Sampling  Per, 


0.1 


DX 


0.001 


DD 


0.00001 


The  Type 
of  The  Data 


Manipulated  to  take  the  step  out 
Rounded  to  seven  significant  digits 


Realized 

System 

Matrix 


0.2220x10  15   -0.1554x10  14   0.2018 


0.1000x10' 


0.3552xl0~14  -0.1104X101 


0.2220x10 


-15 


0.1000x10" 


0.1879x10' 


Realized 
Distribution 
Matrix  T 


0.1000x10 


0.1554x10  14  -0.4440x10  15 


Realized 

Observation 

Matrix 


0.4537 


0.2088x10' 


0.3627x10' 


Eigenvalues 
of  System 
Matrix 


0.9047  +  j  0.3596  x  10 
0.3676  -  j  0.3226  x  10 
0.6067  -  j  0.3617  x  10 


-12 


-12 


-13 


Corresponding 
s -Domain 
Eigenvalues 


-1.001 

•lb.  00  3 

-4.999 
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EXAMPLE  8c 

Transfer 
Function 
of  The 
Sys  tem 

4000 

(s+1) (s+5) (s+10) 

Sampling  Per. 

0.1 

DX 

0.01 

DD 

0.00001 

The  Type 
of  The  Data 

Manipulated  to  take  the  step  out 
Rounded  to  four  significant  digits 

Realized 

System 

Matrix 

0.1332xl0"14    0.1776xl0~14    0.1896 
O.lOOOxlO1      0.0            -0.1080X101 
0.1554xl0~14    O.lOOOxlO1      0.1866X101 

-' 

Realized   

Distribution 
Matrix  T 

O.lOOOxlO1     -O.lllOxlO"14   -0.2220xl0"15 

Realized 

Observation 

Matrix 

0.4537          0.2088X101      0.3628X101 

-14 

-  Eigenvalues 
of  System 
Matrix 

0.8977  -  j  0.1602  x  10 
0.3316  -  j  0.1348  x  10"14 
0.6368  +  j  0.3664  x  10~14 

Corresponding 
s -Domain 
Eigenvalues 

-1.08 
-11.0 
-4.51 
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EXAMPLE  9a 

Transfer 
Function 
of  The 
System 

80000 

(s+40) (s  +4s+400) 

Sampling  Per. 

0.1 

DX 

0.01 

DD 

0.00001 

The  Type 
of  The  Data 

Manipulated  to  take  the  step  out 
Rounded  to  eight  significant  digits 

Realized 
Sys  tern 
Matrix 

0.1942xl0"15   -0.4440xl0~15   0.1227xl0_1 
O.lOOOxlO1     -0.4440xl0~15  -0.6581 
0.0            O.lOOOxlO1    -0.6481 

Realized 
Distribution 
Matrix  T 

O.lOOOxlO1      0.1332xl0~14   0.1332xl0~14 

Realized 

Observation 

Matrix 

0.4355X101      0.3682X101    -0.4645X101 

Eigenvalues 
of  System 
Matrix 

-0.3332  +  j  0.7478 
-0.3332  -  j  0.7478 
0.1831  x  10_1  +  j  0.3330  x  10~15 

Corresponding 
s -Domain 
Eigenvalues 

-1.998  +  j  19.9 
-1.998  -  j  19.9 
-40.0 
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EXAMPLE  9b 

Transfer 
Function 
of  The 
System 

80000 

2 

(s+40) (s  +4s+400) 

Sampling  Per. 

0.1 

DX 

0.01 

DD 

0.00001 

The  Type 
of  The  Data 

Manipulated  to  take  the  step 
Rounded  to  four  significant 

out 
digits 

-0.8881xl0~15   -0.4440xl0"15 

O.lOOOxlO1      0.0 
-O.lllOxlO-14    O.lOOOxlO1 

0.1263xl0_1 
-0.6576 
-0.6474 

Realized 

System 

Matrix 

Realized  , 
Distribution 
Matrix  T 

O.lOOOxlO1      0.6661xl0"15 

0.2220xl0~15 

Realized 

Observation 

Matrix 

0.4356X101      0.3683X101 

-0.4646X101 

Eigenvalues 
of  System 
Matrix 

-0.3331  +  j  0.7478 
-0.3331  -  j  0.7478 
0.1885  x  10"1  +  j  0.5551  x  10"16 

Corresponding 

s-Domain 

Eigenvalues 

-1.998  +  j  19.9 
-1.998  -  j  19.9 
-39.8 
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EXAMPLE  10a 

Transfer 
Function 
of  The 
System 

80000 

(s+40) (s+2+jlOX) (s+2-jl0X) 

Sampling  Per. 

0.1 

DX 

0.00001 

DD 

0.00001 

The  Type 
of  The  Data 

Manipulated  to  take  the  step  out 
Double-Precision 

Realized 
Sys  tern 
Matrix 

-0.2220  x  10~15     0.1499  x  10_1 
0.1000  x  101      -0.8004 

Realized 
Distribution 
Matrix  T 

0.1000  x  101      -0.2220  x  10~15 

Realized 

Observation 

Matrix 

0.2981  x  101      -0.1765  x  101 

Eigenvalues 
of  System 
Matrix 

-0.8187  +  j  0.0 
0.1831  x  10"1  +  j  0.0 

Corresponding 
s -Domain 
Eigenvalues 

-2.0  +  j  10X 
-40.0 
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EXAMPLE  10b 


Transfer 
Func  ti  on 
of  The 
System 


80000 


(s+40) (s+2+jlOX) (s+2-jlOX) 


Sampling  Per, 


0.05 


DX 


0.00001 


DD 


0.00001 


The  Type 
of  The  Data 


Manipulated  to  take  the  step  out 
Double-Precision 


Realized 

System 

Matrix 


-0.8881x10 


-15 


-0.2775x10 


-16 


0.1108 


0.1000x10' 


0.1387x10  16   -0.8187 


-0.4440x10 


-15 


0.1000x10' 


0.1353 


Realized 
Distribution 
Matrix  T 


0.1000x10' 


0.1387x10  16    0.5551xl015 


Realized 

Observation 

Matrix 


Eigenvalues 
of  System 
Matrix 


0.8901 


0.2091x10' 


-0.1333 


0.8471  x  10"7  +  j  0.9048 
0.8471  x  10~7  -  j  0.9048 
0.1353  -  j  0.9048  x  10"14 


Corresponding 
s -Domain 
Eigenvalues 


-2.0  +  j  10X 
-2.0  -  j  10X 
-40.0 
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III.   CONCLUSION 

Ho  * s  Algorithm  provides  a  valuable  tool  for  complete 
identification  of  linear— stationary  systems  which  are  com- 
pletely controllable  and  observable.   The  computations  re- 
quired are  simple  and  can  be  programmed  easily,  in  fact  for 
most  of  the  computations,  available  library  subroutines  of 
a  computer  facility  can  be  used. 

For  discrete  systems  accuracy  of  the  measurement  is 
the  only  limiting  factor.   For  continuous  systems  a  dis- 
crete realization  seems  to  be  more  feasible  as  the  evalua- 
tion of  derivatives  introduces  error.   The  tested  examples 
showed  that  a  five-digits-precision  data  would  give  fairly 
good  results  for  this  approach. 
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DIGITAL    CHMpuTFO    PROGRAMS 


m,a  t  p^    cp  rr, pa M 

pcA.nc;    thf    PAT*    hnjf    ON    EACH    CAPO    w^th 
F^PM^,T(P17.  10)  .     FOP     PI^C^ETF    RE0,P.ES£,NT  ,4T  I ON    FF 
TMr     cySTFM    4T     LPAST    2*M    "AM^LFS     MUST     RC 
MJPDLIFP    pT"cT    BEING    ^A^Lcn    0Nr    FERIOC    AFTFP 
INITIAL    TIME.     php    CONTINUOUS    REPRFSENT AT  ION 
PF     THF     *YrTPM     (EVALUATING    DCP.IV«TTVPS)     «T    LEAST 
o*M  +  l     SAMPLES    MU^T    Be    SU°PLIcO    FIPST    BEING 


S^yPL^n    AT     INITIAL    T!*.Ft     WHFRF 


tuc     CCTIMA  TCD 
FLUC    TISIF. 


IS 


UPPFR  LT^TT  pp  SYSTf.M 


AT    (PAST 
01  M^N*"  ION 


IMDLI 

p  fTM  D  L 

OT  MTM 

G  t  A  A  (  2 

KM  l*c 

r 

MM 
MM=1 

r. 

MP=  1 

c 

Ft 

c 
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PL*G  = 

r 

M 

c 

M=  c 

I  =m*2 

r 

W 

W  =  0.1 

WP  T  Tp 

TF(FL 

c 

PF 

pn    o 
o  F  A  0  ( 

Q 

CPNTT 

c 

PI 
FALL 

GO     TO 

r 

TF 

°1 

DO     12 

PTA0( 

12 

rpNTT 

c 

FI 
H(  1  )  = 

nO    02 

Q2 

H  (  L  L  ) 

r 
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WP  I  TF 

WP  T  TF 
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P.U 
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14 
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CM 

on 
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C.000C1 
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0) 

SAMPLING  PERI0D=»,D16.3,1X, 'SECONDS* , /) 
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SUBROUTINE    SDWFD2 

DURopcE  : 

SUBROUTINE    SDWFD2    FINDS    THE    SUCCESSIVE 
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I  —  jmp  MATRIX  ITS  NORMAL  FORM  WILL  RF  FOUND 

V — ACTU*L    ROW    OTMPM«*TON    OP     A 

M — ACTUAL    COLUMN    OIM^NSTTN    OF    A 

I C — RANK    OF    A 

P--PLFMCNTAPY    ROW    TRANSFORMATION    MATRIX 

q — plCMFNTAPY    COLUMN    TR fl NSFORMAT IPN    MATRIX 

PP  —  I     SUIT^rlf     <^MflLL    POSITIVF    NUMPFR,  OUR  ING 

noFPATION    ANY    NUMBER    WHICH    HAS    AN    ABSOLUTF 

V^LUF    LrSS    THAN    Op    !c     MflOF    FOUAL    TO    7FRO 

M,N,*»,PP,IC     fPF     INPUT«,P     AND    0    ARF    niiTDijTC 


NOTF 


A,c  AND  0  MATPICF^  ARF  PIPST  ARRANGED  IN  THE 
ARP&Y  SUCH  THAT  0  I 5  ON  X0P  AND  P  ON  LEFT  OF 
P    ANP    0    APE    IDENTITY    MATRICFS    AT    THE    REGTNING 


SUR POUT  INF     PHOO(MtN,*  ,DP,P,0,IC) 

IMPLICIT    ocAL*8(A-H,o-7) 

DIMENSION    A(20,20),P(20,2C),0(20,20),X(40,40) 

FIND    THE    01 MFNS IPN    PF    X 
K1=M+N 

INITIALI7F     X    TO    7FRP 

no    1    1=1, ki 

DO     1     J=1,K1 
X(  T  ,J)=0.0 

INTTIALI7E 
K2=N+1 
DO    2    I=K2,K1 


P    CART    OF    X    TO    IDENTITY 


59 


11 


12 
13 
14 


1*5 
16 
17 


IB 
1° 
21 


22 


DP  2  J=l, 
IF( J.NE.f 
X(  !f J)=l. 

CONTINUE 

T  N  T  T 1 1 
K3=M+1 

^o  3  i=i, 

D^  3  J  =  K3 
IF( I.ME.f 
X( T,J)=1. 
CnNTINtJF 

PUT     TH 

00    6    1=1, 

nn    A    J=lt 

X(  (  H-N).( 

K3    IS 

AND    0 

nn    21    J=K 

RIG    IS 

SEAPCH 

Bir,=c.o 

AFTEF 
POW    TP 

IF( ( J-M). 

K4=K2-*3+ 

on    u    i=k 

SEARCH 

MEMORI 

IPfRTG.RF 

bio=pabm 

TPOW=I 
JCOL=J 
CONTINUE 
IPCBICNF 
IF    ALL 
COLUMN 
TO    THC 
00    13    TI= 
TPMP=X( II 
KK1=K1-1 
00     12     JJ= 
X< II,JJ)= 
X(  II  ,K1)  = 
CONTINUE 

Gn   th   5 

IFf  TROW.F 
IF  PIG 
00  15  Jl= 
TrMp  =  Xf  K4 
X( K4, Jl)= 
X( IP0W,J1 
CONTINUE 

DIVIDE 
TFMP=X( K4 

on    it   ji= 

XIK4, Jl)= 

00  1°  I=K 
IF( T.EO.K 
IFf  X(   I  ,J) 

G  =  X< I ,J) 
00     IP     Jl= 
X(  I,J1)=X 

IF(OABS(X 
CONTINUE 
CONTINUE 
CONTINUE 
t  T     thi 

BE     7EP 
TMF    /  , 

TRANSF 
K5=K3+IC 
K6=N+IC 


I-N))  GO  TO  2 

C 

LI7F  0  PART  OF  X  TO  IDENTITY 

N 

,«1 

J-M))  00  TO  3 

0 

E  A  MATRIX  INTO  ITS  PLACE  AT  X  ARRAY 

M 

N 

J+M))=A(I ,J) 

THF    COLUMN    OF    X    WHEPE    FIR^T    COLUMN    ^  F    THE    A 
MATRICES     APE 
3     K 1 

THE    LARGEST    ELFMFNT    IN   COLUMN 

INITIALI ?ED    TO    ZcPO 


IC    COI  IIMNS    OF    A    IS    PUT 

ANSFHPMATIONS    GO    TH    22 

GT.  TO    GO    TP    22 

J 

4,K1 

THE  COLUMN  OP  A  FOR  L 
7E  ROW  AND  COLUMN  OF  L 
.DABSfXU  ,  J)  )  )  GO  TO  I 
X  (  I  ,  J  )  ) 


.O.C)  GO  TO  14 
ELEMENTS  OF  COLUMN  AP 
SHIFT  ALL  COLUMNS  AT 
LFET  AMP  PUT  REMOVEO 

i,Kl 

»J) 


TO  ITS  NORMAL  FORM  PY 


ARREST  FLFMENT  B^LOW  K4 

ARGCST  ELFMENT 

1 


E    ZCRH    PFMPVE    THIS 
RIGHT    OF    THIS    ONE    PLACE 
COLUMN   TO    LAST    COLUMN 


J.KK1 

X(IIt<  JJ+l)  ) 

TEMF 


O.KA)     GO    TO    16 

IS    NOT    IN    APPROPRIATE 
1,K1 
t  II) 

X(IPOWf Jl) 
)=TEMP 

ROW  OF  BIG  BY  BIG 
t  J) 
LrKl 

X(K4,J1)/TEMP 
2»K1 

*)  GO  TO  lo 
.EO.0.0)  GO  TO  19 


ltKl 

f I,J1)-G*X(K4,J1) 

(I, Jl!). IT. DD)  XU,J1)=0.0 


ROW  CHANGE  ROW^ 


S  POINT  ALL  RCWS  OF  TH 
Ot  IF  RANK  EQUAL  TO  CO 
NORMAL  FORM  IS  OBTAINE 
HRMATIONS  FINISH  THE  A 


E  A  BELOW  ICTH  SHOULD 

LUMN  DIMENSION  OF 

0,  OTHERWISE  BY  COLUMN 

LGORITHM 


60 


31 
32 

33 

*■  1 

23 


24 


r 
r 
c 
c 
c 

r 

r 
C 

c 
c 
c 

r 

r 
C 

c 

r 
c 
c 


U=K3 

ip(  ir  . 

OP  33 
nn  32 
IF(  Xf  L 

r,  =  x(L, 

on    31 


Fn.fv)    GO    TO    41 
I  =  *  ?  ,  K  6 
|  =  KS,K1 
,J) .FO.O.O) 
J) 
1=1, Kl 


OP    TO    32 


ITS    ARRAY 


xc T ,  J)  =  x(  f  ,.n-o*x(  I  ,i  l ) 

!F(0*rc(  Xf  T  ,  J)  )  . LT.PP)     X(I  ,  J)  =0.0 

CONTINUE 

CONTINIJF 

L  1  =  1  1«-1 

CONTINUE 

PUT  P  to 
nn  ?3  1=1,  m 
OP  23  J  =  l,v 
°f  t  ,  J)  =  Xf  II  +  M  t  J) 

PUT    o    TO    ITC     ARRAY 
OP    2^     1=1, N 
00    24    .1=1, N 

Of  I,J)=Xt  I  ,f  J+*)  ) 

DPT||PM 

FNP 


SUBROUTINE    mult 

PURPOSF  : 

SUBROUTINE     Mill  " 

corner  pf  c=ab 


FINOS    THF    (IC.IM)    NORTHWEST 
MATRIX     mijltIPL  IC^TTPN 


nFSrPIPTITN    PF     PAPAMETCRS: 
£  —  MULTIPLICAND    MATRIX 
B — WULTIPLTCP    V4TPIX 
f — RESULTANT    MATRIX 

1 — ACTUAL    rri"MM    OIMcnjct^N    HP    THC     A 
If — crw    DIMENSION    OF    OEMREP    NORTHWE*T    CORNFR 
!M — COLUMN    OIMENSIPN    OF    DESIRcO    N.W.    CO»NFR 
A,R,L,IC,IM    APE    INPUTS, C    IS    OUTPUT 

SUBROUTINE    MtLTf  L,A,B  ,IC,IM,C) 

TMCLITIT    PEftL*F( A-H, D-7) 

DIMENSION  A f 20,20)  ,B f 20  ,20 )  ,C ( 20  ,20  ) 

OP  3  1=1, IC 

np  2  J=1,IM 

$um=o. c 

OP  1  K=1,L 

cUM=C|)y  +  A  (  i  tk)  *P(  K,  ,1) 

Cf I,J)=SUM 

CONTINUE 

RETURN 


61 


BIBLIOGRAPHY 


1.  B.  L.  Ho  and  R.  E.  Kalman,  "Effective  Construction  of 
Linear  State-Variable  Models  from  Input/Output  Func- 
tions", Sonderdruck  aus  "Regelungstechnik",  14.  Jahrg. , 
Heft  12,  1966,  Seite  545-548. 

2.  R.  E.  Kalman,  "Mathematical  Description  of  Linear 
Dynamical  Systems",  J. S.I. A.M.  Control,  Ser.  A,  Vol.1, 
No.  2,  1963. 

3.  Elmer  G.  Gilbert,  "Controllability  and  Observability 
in  Multivariable  Control  Systems",  J. S.I. A.M.  Control, 
Ser.  A,  Vol.  2,  No.  1,  1963. 

4.  Paul  Zeiger,  "Lecture  Notes  on  Algebraic  Theory  of 
Linear  Systems". 

5.  Robert  C.  K.  Lee,  "Optimal  Estimation,  Identification, 
and  Control",  Research  Monograph  No.  28,  The  M.I.T. 
Press. 

6.  R.  E.  Kalman,  "On  the  General  Theory  of  Control  Sys- 
tems", Proc.  1st  International  Congress  on  Automatic 
Control,  Moscow,  1960. 


62 


INITIAL  DISTRIBUTION  LIST 


No.  Copies 


1.  Defense  Documentation  Center  20 
Cameron  Station 

Alexandria,  Virginia  22314 

2.  Library,  Code  0212  2 
Naval  Postgraduate  School 

Monterey,  California  93940 

3.  R.  E.  Kalman  1 
Department  of  Operation  Research 

Stanford  University 
Stanford,  California  94305 

4.  Asst.  Professor  D.  E.  Kirk  1 
Department  of  Electrical  Engineering 

Naval  Postgraduate  School 
Monterey,  California  93940 

5.  Asst.  Professor  J.  S.  Demetry  1 
Department  of  Electrical  Engineering 

Naval  Postgraduate  School 
Monterey,  California  93940 

6.  Assoc.  Professor  H.  A.  Titus  1 
Department  of  Electrical  Engineering 

Naval  Postgraduate  School 
Monterey,  California  93940 

7.  Professor  G.  J.  Thaler  10 
Department  of  Electrical  Engineering 

Naval  Postgraduate  School 
Monterey,  California  93940 

8.  Dr.  R.  C.  Dorf,  Dean  1 
School  of  Engineering 

Ohio  University 
Athens,  Ohio  45701 

9.  LT  I.  Eldem  4 
Giinaydin  Mah. 

Istiklal  Cad  Ar-is   Ap.  Kat.3 
Bandirma,  TURKEY 

10.   Deniz  Kuvvetleri  Komutanligi  1 

Ankara,  TURKEY 


63 


No.  Copies 


11.  Deniz  Harb  Okulu  Komutanligi  1 
Heybeliada  -  Istanbul  -  TURKEY 

12.  Middle-East  Technical  University  1 
Ankara,  TURKEY 

13.  Istanbul  Teknik  Universitesi  1 
Taskisla  -  Istanbul,  TURKEY 

14.  Karadeniz  Teknik  Universitesi  1 
Trabzon,  TURKEY 


64 


Unclassified 


Security  Classification 


DOCUMENT  CONTROL  DATA  -R&D 

,  Security  clas  sitication  ol  title,   body  ol  abstract  and  indexing  annotation  niu.vl  he  entered  when   the  overall  report  Is   classified 


I     originating    ACTIVITY  (Corporate  author) 

Naval  Postgraduate  School 
Monterey,  California  93940 


2«.   REPORT    SECURITY    CLASSIFICATION 


Unclassified 


2b.  GROUP 


3   REPORT  TITLE 


Identification  of  Plant  Dynamics  Using  Ho's  Algorithm 


4.    DESCRIPTIVE  NOTES  (Type  ol  report  and.inclusive  dates) 

Master's  Thesis;  June  1969 


5     authorisi  fF/rsf  name,  middle  initial,  last  name) 


Ilker    ELDEM 


6      REPOR  T    DATE 


June    1969 


la.    TOTAL    NO.    OF   PAGES 

66 


7b.    NO.    OF    REFS 

6 


3a.    CONTRACT    OR    GRANT    NO. 


b.    PROJEC  T    NO. 


9a.    ORIGINATOR'S    REPORT    NUMBER(S) 


9b.   OTHER   REPORT  NOISI  (Any  other  numbers  that  may  be  assigned 
this  report) 


10.    DISTRIBUTION    STATEMENT 


Distribution  of  this  document  is  unlimited. 


11.    SUPPLEMENTARY    NOTES 


12.    SPONSORING   MILI  TARY    ACTIVITY 


Naval  Postgraduate  School 
Monterey,  California  93940 


13.  ABSTRACT 


Ho's  Algorithm  is  reviewed  and  demonstrated  with  analytic  exam- 
ples.  A  digital  computer  program  is  developed  to  implement  the  algo- 
rithm for  single  input,  single  output  systems  and  used  to  identify 
linear  continuous  and  stationary  systems  which  are  driven  with  a  unit 
step  as  the  test  input.   Discrete  realization  of  the  continuous  sys- 
tems is  obtained  using  the  measured  output  samples,  to  a  step  input, 
directly  in  the  algorithm. 


DD 


F0B»  1473 

1    NOV   »B   It    /    W 

S/N    0101 -807-681 1 


(PAGE     1) 


Unclassified 


65 


Security  Classification 


A- 3 1408 


Unclassified 


Security  Classification 


KEY     WORDS 


LINK     A 


LINK     B 


Identification 
Algorithm 
Linear  Systems 


DD  ,Frr..1473  'back 


/'J    0101-907-6821 


Unclassified 


Security  Classification 


66 


thesE312 

Identification  of  plant  dynamics  using  H 


II  MUM  II 
3  2768  000  98505  5 

DUDLEY  KNOX  LIBRARY 


